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qv ABSTRACT 
On 

We explore the dynamical restrictions on the structure of dark matter halos through 
Q-r a study of cosmological self-similar gravitational collapse solutions. A fluid approach 

■ to the collisionless dynamics of dark matter is developed and the resulting closed 

, set of moment equations are solved numerically including the effect of halo velocity 

dispersions (both radial and tangential), for a range of spherically averaged initial 
density profiles. Our results highlight the importance of tangential velocity dispersions 
to obtain density profiles shallower than 1/r 2 in the core regions, and for retaining 
qq ' a memory of the initial density profile, in self-similar collapse. For an isotropic core 

CN ■ velocity dispersion only a partial memory of the initial density profile is retained. If 

■ tangential velocity dispersions in the core are constrained to be less than the radial 
ON . dispersion, a cuspy core density profile shallower than 1/r cannot obtain, in self-similar 

collapse. 

2 ' Subject headings: Cosmology: dark matter, Large-scale structure of Universe; Galaxies: 

c/2 . Formation, Halos, clusters 



1. Introduction 

In hierarchical clustering theories of structure formation, like the cold dark matter (CDM) 
models, small mass clumps of dark matter form first and gather into larger and larger masses 
subsequently. The structure of these dark matter " halos" , is likely to be related to how the halos 
formed, the initial spectrum of the density fluctuations and to the underlying cosmology. Several 
properties of galactic and cluster halos can be well constrained by observations. So, if the matter 
distribution in dark halos are fossils which do depend on some of the properties of structure 
formation models, like their initial power spectrum, one would have a useful observational 
handle on these properties. It is therefore necessary to understand what determines the matter 
distribution (or density profiles) of dark matter halos ab initio. This forms the motivation of the 
present paper and the companion paper by Subramanian, Cen and Ostriker (SC099, 1999). 

Further, Navarro, Frenk and White (NFW) (1995, 96, 97) have proposed from their N-body 
simulations, that dark matter halos in hierarchical clustering scenarios develop a universal density 
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profile, regardless of the scenario for structure formation or cosmology. The NFW profile has an 
inner cuspy form with the density p oc r _1 and an outer envelope of the form p oc r -3 . Several 
investigators have found that the NFW profile provides a moderately good fit to numerical 
simulations ( Cole and Lacey 1996, Tormen, Bouchet & White 1997, Huss, Jain & Steinmetz 1997, 
1999, Thomas et al, 1998). Recently, though, high resolution simulations of cluster formation in 
a CDM model, by Moore et al. (1998), yielded a core density profile p(r) oc r~ 1A , shallower than 
r -2 , but steeper than the r _1 form preferred by NFW, consistent with the earlier high resolution 
work of Xu (1995). (For smaller mass halos, Kravtsov et al. (1998) find the core density profile to 
be shallower than the NFW form). It is important to understand these results as well on general 
theoretical grounds. 

In a companion paper SC099, we explore the possibility that a nested sequence of undigested 
cores in the center of a halo, which have survived the inhomogeneous collapse to form larger and 
larger objects, determine halo structure in the inner regions. For a flat universe with a power 
spectrum of density fluctuations P(k) oc k n , scaling arguments then suggest that the core density 
profile scales as, p oc r~ a with a = a n = (9 + 3n)/(5 + n). But whether such a scaling law indeed 
obtains depends on the detailed dynamics. 

Similarity solutions often provide a tractable, semi-analytic route to study time dependent 
dynamics in complicated physical systems. Fillmore and Goldreich (FG, 1984) and Bertschinger 
(B85, 1985) derived such solutions for describing the purely radial collapse of cold, collisionless 
matter in a perturbed Einstein-de Sitter universe. These solutions need to be generalised to 
incorporate tangential velocity dispersions, which as we see below, turn out to be crucial to 
understand density profiles shallower than 1/r 2 . Some general, analytical aspects of the similarity 
solutions incorporating tangential velocity dispersions are outlined in the companion paper SC099. 
In the present paper, we consider these self-similar collapse solutions in greater detail, by deriving 
and solving numerically the scaled moment equations for such a collapse, including the effect of 
velocity dispersions. 

In the next section we formulate the self-similar collapse problem and introduce a fluid 
approach for its solution, recapitulating the corresponding discussion in SC099. In Section 3 
we derive scaled moment equations describing the collapse and discuss their numerical solution. 
Specific numerical examples of self-similar collapse are given in Section 4. These solutions include 
the effect of halo velocity dispersions (both radial and tangential), and consider a range of 
spherically averaged initial density profiles. The final section discusses the results and presents 
our conclusions. 

2. The self-similar solution 

We summarize in this section, some of the properties of the similarity solution, that can be 
derived by analytic arguments. Although much of this section is mostly a recapitulation of Section 
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3 of SC099, we include it here to make the present paper as self-contained as possible, and also to 
set the framework for the detailed numerical work which follows. 

Consider the collapse of a single spherically symmetric density perturbation, in a flat 
background universe. Assume the initial density to be a power law in radius. We expect to 
describe the dynamics through a self-similar solution. FG and B85 looked at the self-similar 
evolution by following the self-similar particle trajectory. We adopt a different approach, by 
examining directly the evolution of the phase space density. During the course of this work we 
have learned that several authors (Padmanabhan 1994, unpublished notes; Padmanabhan 1996a, 
Chieze, Teyssier and Alimi 1997; Henriksen and Widrow 1997) have also adopted this approach to 
the purely radial self-similar collapse of FG and B85. We will emphasise and incorporate here also 
an additional aspect, the distinctive role of non-radial motions (velocity dispersions) in self-similar 
collapse. 

The evolution of dark matter phase space density /(r,v,i) is governed by the Vlasov 
Equation, 

df df df 

where r and v = r are the proper co-ordinate and velocity of the particles respectively. Also the 
acceleration a = v = — V<&, with 

V 2 $ = AttG P = 4ttG J fd 3 w (2) 

By direct substitution, it is easy to verify that these equations admit self similar solutions of the 
form 

/(r, v, t) = k 2 k^t-^F(^, ^); p = q + 1 (3) 

where k\,k 2 are constants which we will fix to convenient values below. We have used proper 
co-ordinates here since the final equilibrium halo is most simply described in these co-ordinates. 
(The same solution in co- moving co-ordinates is given in Padmanabhan (1996a)). Defining a new 
set of co-ordinates y = r/(kit p ), w = v/(kit q ) and a scaled potential x = k~^ 2 t~ 2q <&, the scaled 
phase space density F satisfies 

dF dF dF „ dF 

-(q + 2 P )F - P y.— - ,w.— + w.— - V yX .— = 0; (4) 

V 2 yX = 4nGk 2 J Fd 3 w. (5) 

Consider the evolution of a spherically symmetric density perturbation, in a flat universe 
whose scale factor a(t) oc t 2 ^ 3 . For self similar evolution, the density is given by 

p(r, t)= J fd\ == k 2 r 2 J F(y, w)d 3 w = k 2 r 2 ^{y) (6) 

where we have defined r = |r|, y = |y| and used the relation p = q + 1. For the flat 
universe, the background matter density evolves as pi,(t) = l/(6irGt 2 ). So the density contrast 
p(r,t)/pb(t) =ip(y), where we take k 2 = 1/(QttG). 
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2.1. Linear and non-linear limits 

Let the initial excess density contrast averaged over a sphere of co-moving radius 
x = r/a(t) oc rt~ 2 l 3 be a power law 5(x,U) oc x~ 3e . Since p/pb is a function of y alone, the 5{x,t) 
will also be a function only of y. Note that, in the linear regime, the excess density contrast 
averaged over a co-moving sphere, grows as the scale factor a{t). So one can write for the linear 
evolution of the spherical perturbation 

5(r, t) = 5 x- 3 H 2 / 3 oc 5 Q r- 3 H 2 ' 3+2 * oc 5 o2 T 3e r 3ep+2 / 3+2e , (7) 

where we have substituted r oc yt p . This can be a function of y alone, for a range of t in the linear 
regime iff — 3ep + 2/3 + 2e = 0, which gives 

2 + 6e 

P = — (8) 

We see that once the initial density profile is specified the exponents p, q of the self similar solution 
are completely determined. 

Consider now what happens in the non- linear limit. The zeroth moment of the Vlasov 
equation gives 

^ + V r .(pv) = (9) 

Here v =< v > is the mean velocity. (Henceforth both <> or a bar over a variable denotes a 
normalised moment over /). In regions which have had a large amount of shell crossings, it seems 
plausible to demand that the halo particles have settled to nearly zero average infall velocity, that 
is v r = 0. From (^) , we then have (dp/dt) = 0, and therefore, in the non-linear regime, 

p(r, t) = Q(r) = Q{yt p ) = g^^(y)- (10) 

This functional equation has only power law solution, because of the power law dependences on t. 
Substituting Q(r) = qor~ a into Eq. (|l0|), and using r oc yt p , we obtain y- a t- pa oc t~ 2 D(y). This 
can only be satisfied for range of t in the non-linear regime provided pa = 2. So, for an initial 
density profile with a power law slope 3e, the power law slope of the density in the non-linear 
regime is given by, 

*-* = -*-. (id 

p 3e + l v ; 

This result has been obtained by following the self similar particle trajectory, by B85 (for e = 1), 
and FG for 2/3 < e < 1. We see that it can be simply obtained by just combining the self-similar 
solution / and the static core condition. (Obtaining the B85/FG result in this way has been 
independently noted by Padmanabhan (private communication, unpublished notes 1994)). 

What should we choose for the value of e? For a power law P(k) oc k n , the fractional density 
contrast averaged over a co-moving sphere of radius x, is distributed as a Gaussian, with a variance 
oc x -( 3 + n )/ 2 . This suggests a "typical" spherically averaged initial density law for a halo collapsing 
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around a randomly placed point of the form 5(x,ti) oc x ( 3 + n )/ 2 5 or 3e = (3 + n)/2. Suppose we 
use this value of e for the initial density profile of a halo. Then the halo density in the staic core 



regions will be p(r,t) oc r a , where, substituting 3e = (3 + n)/2 in Eq. ( 11 



9 + 3n 

a = a n = — (12) 

5 + n 

Remarkably, this is the same form that scaling laws suggest, for the core of a collapsed halo, 
assuming that the cores of sequence of sub-halos are left undigested, during the formation of the 
bigger halo (see SC099). (In a paper which appeared during the course of this work, Syer & 
White (1998) motivate the same form, in the case when bigger halos form by purely merger of 
smaller halos). Note that for n < 1 the density law given by ( |l2|) is shallower than 1/r 2 . 

FG also showed that a power law slope shallower than 1/r 2 , cannot obtain for purely radial 
collapse, And that while the above form for a should obtain for 2/3 < e < 1, for e < 2/3, one goes 
to the limiting value a = 2. However, this is only true for purely radial trajectories (cf. White 
and Zaritsky 1992; Sikvie, Tkachev and Wang 1997). We see below, by considering the higher 
moments of the Vlasov equation, that a < 2 can only obtain if the system has non-radial velocity 
dispersions. 



2.2. Jeans and Energy equations 

Suppose we multiply the Vlasov equation by the components of v and integrate over all v. 
Assume there is no mean rotation to the halo, that is v$ = and = 0. Then we get 

^ + ^ + >?-^) + ~ = (13) 

») = »J (") 
Here M(r) is the mass contained in a sphere of radius r. 

Let us consider again a static core with v r = 0. The Jeans equation gives two equations 
for the three unknown velocity dispersions, even for a static core. To see if one can close the 
system SC099 considered the second moments of the Vlasov equation (the energy equations). 
However these will involve the third moments, or the peculiar velocity skewness. Some form of 
closure hypothesis is needed in a fluid treatment of the Vlasov equation. For this we proceed 
as follows: One can firstly assume that initially the tangential velocities have zero skewness. 
Then in purely spherically symmetric evolution they would not develop any skewness, that is 
Vq = =< vqv^ >= for all times. Also if the initial velocity ellipsoid had one of its principle axis 
pointing radially, we do not expect this axis to become misaligned in purely spherical evolution. 
This means we can assume < v r Vg >= v r Vg. Under these assumptions, and taking the static core 
condition v r = 0, we get, (d(pvg) / dt) = or pvg = K(r) independent of t. For the self-similar 
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solution we then have 

P v e 



"| = K{r) = K{yt*) = k 2 k 2 t A ^ J w 2 e F(y, w)d 3 w (15) 



Once again substituting a power law solution K(r) = Kor s , to this functional equation, we get 
the constraint from matching power of t on both sides, ps = Aq — 2p. Using p = q + 1, we then get 
s = 2 — 4/p = 2 — 2a, and so 

pvj = K r 2 ~ 2a (16) 



Integrating the radial momentum equation using Eq. (|l^) , (14), (16) and using p = q^r 
we have 



v? = r 2 - a 



K Q 4irGqo 



1 



(2-a) 



{2-a)q 2(2-a)(3-a) 
v 2 (r) 



,2/^ GM(r) 



2r 



(17) 



Several important points are to be noted from the above equation 1 . A crucial one is that, when 
a < 2, the RHS of Eq. (|17|) can remain positive, provided one has a non zero tangential velocity 
dispersions. If one has a purely spherically symmetric collapse and zero tangential velocities, then 
the density law cannot become shallower than a = 2 and maintain a static core with v r = 0. This 
agrees with FG. Infact for any a < 2, one needs tangential velocity dispersions to be at least 
as large as GM/2r, comparable to the gravitational potential energy per unit mass. Further, 
one can see that to obtain static cores with a < 1, the required tangential dispersions have 
to be necessarily larger than the radial velocity dispersions. Also note that for a < 2, all the 
components of velocity dispersions decrease with decreasing radius, as suggested by the simple 
scaling arguments of SC099. 

For a static core v% should be independent of t. However suppose we look at the the energy 
equation for the radial velocity dispersion, 

d(ov 2 ) 1 d(pr 2 v 3 ) 2p < v r {vl + vV> > 

uyfw^ + i u { pr u r) _ _y_ ^_ + 2v ~ rpGM/r 2 = Q lg) 

at r l or r 

This shows that, even when v r = 0, a time independent radial velocity dispersion can only obtain 
if the radial velocity skewness < (v r — v r ) 3 > is also zero. In the core regions where large amounts 
of shell crossing has occurred, one can assume that a quasi "equilibrium" state obtains, whereby 
all odd moments of the distribution function, over (v — v), may be neglected. Such a treatment 
will correspond to considering a fluid like limit to the Vlasov equation. 

However, the radial skewness will become important near the radius, where infalling matter 
meets the outermost re-expanding shell of matter. This region will appear like a shock front in the 



A constant of integration could have arisen in integrating the Jeans equation over radius; however this is excluded 
as, for a static core, arguments similar to deriving jl5| ) and (|lrj), give pv^ oc r 2 ~ 2a , for the self-similar solution. 
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fluid limit. A possible treatment of the full problem in the fluid approach to the Vlasov equation 
then suggests itself. This is to take the radial skewness to be zero both inside and outside a 
"shock or caustic" radius, whose location is to be determined as an eigenvalue, so as to match the 
inner core solution that we determine in this section with an outer spherical infall solution. One 
has to also match various quantities across this "shock", using jump conditions, derived from the 
equations themselves. To do this requires numerical solution of the self consistent set of moment 
equations derived from the scaled Vlasov equation (the main focus of this paper), to which we 
now turn. 



3. Numerical solution of moment equations for self similar collapse 



3.1. Moment equations 

We write the scaled Vlasov equation (0) in spherical co-ordinates and take moments. Let 



us define V 



w r 



II =< (w r — w r ) 2 > and £ 



Wn 



We also set the tangential velocity 



skewness to zero. As explained above, we take the radial skewness to be zero both inside and 
outside a "shock or caustic" radius. The shock location, say y = y s in scaled co-ordinates, will 
be determined as an eigenvalue, to the complete problem. So we set < (w r — w r ) 3 >= in the 
regions of interest, y < y s and y > y s . The resulting moment equations can be further simplified 
with a little algebra and then can be written in the following more transparent form: 



y 2 dy 



y 2 ipV 



2ip -py 



d^ 
dy 



(p - l)ipV + i/)(V 



(V -py)- 




o 



+ 



Mtp 



y 



6iry 2 



d 



dy L 

dM 
dy 



\(ipy 2 ) 3 

+ (4p - 2)Sy 2 



2p- 2 



(19) 

(20) 
(21) 
(22) 
(23) 



These equations have obvious meaning: Eq. (|19|) is the continuity equation for the scaled 
density, (|20| ) the scaled Euler equation and ( |2~T|) the scaled energy equation. Angular momentum 
conservation is reflected in Eq. (|22|) for S. It should be noted that the energy equation reflects the 



more general conservation P/A 3 along fluid trajectories, where A = pr 2 is an effective linear radial 
density and P = A < (v T — v r ) 2 > is the effective radial pressure. Infact our system behaves like a 
mono dimensional gas with an effective adiabatic index 7 = 3, provided one takes the density to be 
the linear radial density A, and defines the pressure P as above. 
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Both radial and tangential velocity dispersions are likely to be generated during the 
inhomogeneous collapse to form the halo. So it is natural take the initial, pre-collapse, velocity 
dispersions to be small. In the problem we are treating, of purely spherical collapse, the radial 
velocity dispersions will automatically be generated when spherically collapsing shells start to 
cross, that is where the radial skewness is important. In our fluid approach we will be replacing 
this region where radial skewness is important by a shock front. On the other hand, note that 
there are no source terms in the moment equations for the tangential velocity dispersions. Indeed 
in a purely spherically symmetric problem tangential velocities have to be necessarily introduced 
in an ad-hoc fashion. They are only non zero if present in the initial conditions. 

White and Zaritsky (1992) introduced tangential velocities into their solutions by invoking 
a fictitious tangential force which act on particles in each spherical shell, until the shell turns 
around. Since in a general inhomogeneous collapse, one expects all the components of the velocity 
dispersion to be generated together, the shock gives an alternate natural location to introduce a 
tangential velocity dispersion, as well. We will do this here, for most of the numerical examples. 
Further in a non spherical, asymmetric collapse, the random velocities induced at the shock 
would be in general non-radial (cf. Ryden 1993), and the spherical models with both tangential 
and radial velocity dispersions introduced at the shock, may represent this in a rough way. (For 
comparison, we will also present a few examples in the next section, with the tangential velocity 
dispersion introduced at the turn around radius.) Studies of halo formation using cosmological 
N-body simulations, which are discussed in SC099, can redress this deficiency of the spherical 
treatment. 

The evolution of the region before shell crossings is determined by the spherical infall solution. 
At some initial time ti, let the excess density contrast averaged over a sphere of proper initial 
radius r{ be <5j(rj) = <5o(rj/ro)~ 3e . Then the shell initially at r, will turnaround and collapse when 

— — o /o 

it has expanded to a radius rj/<5j(rj) at a time t = (3^/4)^/^ (r^). The radius of a shell turning 
around at any time t is given by r t (t) = ro t (t/tot) p where p = (2 + 6e)/9e is as in Eq. (||) . Also 
r ot = (^o/^o) and iot = (3^/4)^ /5 > are the turn-around radius and time of the shell initially at 
vq. Since y = r/(kit p ), a natural way of fixing the constant k\ is by taking k\t p = rt{t). We will 
do this in what follows. Then turn around occurs at the scaled co-ordinate y = 1. 

A straightforward application of the spherical model (cf. Peebles 1980,, Padmanabhan and 
Subramanian 1992, Kumar, Padmanabhan, Subramanian 1995) then gives the solution of the 
moment equations, when II = S = 0, in the region y > y s . Expressed in a parametric form we 
have for y > y s , 



r 




V(y) = y 




II 



r t (t) 2(0-sin0)P 




(24) 



This goes over to the standard growing mode solution in the linear limit as y — > oo. 
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3.2. Matching and boundary conditions 

The equations ( p4| ) evaluated at y = y s gives the pre-shock boundary conditions to the 
moment equations. To match the spherical infall solution to the core solution for y < y s determined 
in the previous section, we have to specify the jump conditions across the shock at y = y s . These 
conditions can be derived again in a straightforward manner from the moment equations. Suppose 
we denote the pre-shock values with subscript 1 and post shock values of all quantities with a 
subscript 2. Also we wish to consider the case when the pre-shock IT = 0. Then the scaled jump 
conditions are given by 

ih = 2ih, V 2 =py s + ~(V 1 -py s ), n 2 = (Fl ~ PVs)2 , M 2 = M 1 (25) 

Infact the jump conditions corresponds to taking 7 = 3 in the usual fluid Rankine-Huginot jump 
relations. These together with a non zero, arbitrary £2, gives the starting values for the numerical 
integration of the scaled moment equations fllS| ) - (22) inward from the shock location y = y s . The 



eigenvalue y s is determined by requiring the solutions to satisfy the inner boundary conditions 

V = M = 0, y = (26) 

To ensure the vanishing of the mass at y = 0, we have in fact integrated the scaled continuity 
equation and expressed the scaled mass interms of the density and velocities. We have using (19) 



and (23) 



M( V) _ (27) 

The scaled density for all a and the scaled dispersions for a > 2, are expected to be singular at 
the origin for the shocked infall solutions. So we scale out the expected asymptotic behaviour, at 
y — > 0, before numerical integration. If we are to obtain a nearly static core, we expect V — > and 
dV/dy — > as y — > 0. In this case, an analysis of the moment equation shows (see also section 2), 

1>(v) = y~ a i>(y), EG/) = y^iv), n( y ) = y 2 ~ a ii{y) (28) 

where tp(y), £(y) and n(y) are expected to tend towards a constant value as y — > 0. The exact 
asymptotic dependence of V(y) of course has to be determined by the numerical solution. 

The moment equations ( |i~9| ) - ( |22| ) are numerically integrated, after eliminating the scaled 
mass using ( |27| ) and transforming to the dependent variables defined in Eq. . We adapted 

a NAG library routine which integrates the differential equations using a Runge-Kutta-Merson 
method, and solves the boundary value problem with Newton iteration in a shooting and matching 
technique. For a given e, and a sufficiently large £2 (when a < 2), a unique value of y s is found 
to satisfy the inner boundary conditions of (^). The moment equations lead to two conservation 
laws which can be used to provide a check on the numerical integration. These can be derived by 
using Eq. ( p7|) and the moment equations. We have 

- yy! = const, k = (29) 

M K (y) 3- a v ; 
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representing angular momentum conservation, where M{y) = y s a M(y). And 

Sj»^M = const, „=§^ (30) 

representing energy conservation. At each point these two integrals of motion were checked and 
had relative errors less than 10 -8 — 10 -9 , as in B85. A possible additional constraint on a solution 
is the asymptotic condition given by Eq. (|l7|), for an almost static core. In terms of the scaled 
variables, static core solutions satisfy the constraint 

-0(0) 



n(o) 1 



(2 -a) 



S(0) 



3(3 - a) 



(31) 



The equations of course cannot be integrated in practice upto y = 0, but only upto some small 
V = Um which is generally ~ 10~ 2 ' 5 — 10 -4 . So this constraint can be checked at this minimum y. 
Also, in practice, V(y m ) is not expected to be identically zero (unlike V(0)), and so one has to set 
it to a very small but non-zero value to obtain converged solutions. In general for the solutions 
obtained here, V(y m ) = V m ~ —1.0 x 10~ 6 to —1.0 x 10~ 3 and V m /y m « 1. And the constraint 
Til) is satisfied at a few percent level. Let us now discuss particular numerical examples. 



4. Numerical examples 

4.1. Collapse onto a localised, overdense perturbation, e = 1 case 

First we look at self-similar spherical secondary infall onto an initially localised, overdense 
perturbation, by adopting e = 1 and £2 = 0. This problem was solved by B85 and FG by 
examining the self similar particle trajectory. The parameters for the numerical solution obtained 
here with the fluid approach are summarized in Table 1. We give there the assumed parameters, 
the derived eigenvalue y s and the value of the scaled dependent variables at y s and at the 
minimum y = y m . We find the eigenvalue y s = 0.4628 for the above parameters. B85 solving the 
problem by looking at particle trajectories got the location of the outermost caustic as y s = 0.364. 
The difference between the location of the shock as determined by this work and B85 could be 
because we have replaced a smooth transition region for the collisionless fluid, where velocity 
skewness is important, by a discontinuous shock. B85 found that the scaled density could be fitted 
asymptotically by a form ip(y) ~ 2.79y _9//4 when they adopted a minimum y = y m ~ 0.02 for 
the particle trajectory. We can integrate our equations and get converged solutions satisfying the 
boundary conditions upto y m ~ 2 x 10~ 4 . We find that ip(y) = t/jy~ 9 / 4 ~ 3.1y~ 9 / 4 at y m , while at 
y ~ 2 x 10 -2 we find ip ~ 2.5. These numbers bracket the asymptotic value of tp ~ 2.79 obtained 
in B85. So there is reasonable agreement between our work and B85, given the differences in the 
value of y m and the very different approaches. 

In Figure 1 we plot V(y), log{ip{y)) } log(Tl(y)), against log{y) for this solution with e = 1, 
T,2 = 0. We can define the scaled rotational velocity U(y) = [(GM(r)/r)/(r^/t) 2 ] 1//2 . In Figure 
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1 we also show a plot of log(U 2 ) versus log(y). We see that the velocity, V(y), smoothly tends 
to zero as y — > 0. V m for this case was —4.0 x 10~ 5 . To compare the asymptotic dependence 
of the scaled density, with that predicted above for a static core, we also show in the log{ip) 
vs log{y) plot, density laws ip(y) oc y~ a , with a = 9e/(l + 3e) predicted in section 2, (dashed 
line) and ip{y) oc y~ 2 (dot-dashed-dot line). These are normalized to agree with ip(y) at the 
minimum y shown in the figure. We see from figure 2 that as y — > 0, the density does go over to 
ip(y) oc y~ a oc y~ 9//4 , as expected for a static core, with e = 1. Overall, we recover the results of 
B85 and FG reasonably well with our fluid approach to the problem. 



4.2. e < 2/3 cases and the importance of tangential dispersions 

We then considered solutions for initial density profiles shallower than r~ 2 , or e < 2/3. For 
such shallow density profiles, if the collapse were purely radial, FG showed that the final density 
profile approaches a 1/r 2 form. We find, as expected, that the nature of the solutions in this case, 
depends on the ratio of tangential to radial velocity dispersions. We illustrate this by considering 
two values of this ratio, which bracket the expected behaviour. 

In Figure 2 we show the solution for the case e = 0.4, £2 = ^2Ul~ a = 0.94. The detailed 
solution for this case is given in Table 2. For this solution, the value of y s = 0.4955. We show 
both log(Ii{y) (solid line) and log(Ti{y) (dashed line) in the same plot, so that they can be easily 
compared. From the figure or the table one sees that tangential velocity dispersions for this 
solutions are everywhere larger than the radial dispersions, by a factor ~ 1.3. (Or (£/n) 1//2 ~ 1.3). 
For e = 0.4, and a static core, we expect the scaled variables, to have the asymptotic behaviour 



given in Eq. fl28|) , with a = 18/11. We see from comparing the solid and dashed lined, in 
the ip(y) — y plot of Figure 2, that the density rises as ifi oc y~ a oc y~ 18 / n to a very good 
approximation, throughout the core. Also the velocity dispersions and rotation velocities decrease 
with decreasing radius, as the analytic theory of Section 2 (or Eq.(f28|)) predicts. Indeed, from 
Table 2, we see that all the variables tp(y), S(y) and n(y) tend to constant values as y — > to an 
excellent approximation. This case illustrates that it is possible to obtain solutions for e < 2/3, 
which have a = 9e/(l + 3e) < 2, provided the tangential velocity dispersions are large enough. 

To illustrate the effect of decreasing tangential velocity dispersions, we show in Figure 3, 
the properties of a solution with e = 0.4, £2 = 0.65. The parameters for this solution, are given 
in Table 1. The location of the shock is at y s = 0.3797. We could get converged solutions with 
Vm = 4.6 x 10~ 3 , and with the velocity V m = —1.75 x 10~ 3 . The core regions are nearly static 
but not completely so. But the constraint given by Eq. (|3~l|) is satisfied at the 2% level. For 
this case the radial velocity dispersions are everywhere larger than the tangential dispersions, by 
a factor ~ 1.15 as y — > 0. One sees a large difference between this solution (Figure 3) and the 
one obtained for larger tangential velocity dispersion (Figure 2). First we see that when radial 
dispersion dominates, the density profile is closer to the rp oc y~ 2 form than the tp oc y~ a form, 
although neither provides a good fit. Second the velocity dispersions are reasonably constant with 



-12- 



radius as y — » limit, instead of decreasing with decreasing radius. The rotation velocity is also 
flatter. 

We also considered smaller values of e. Figure 4 gives a solution with e = 1/6, (corresponding 
to a = 1 or n = —2 in a n ), and S2 = 2.2 and some parameters for this solution are given in 
Table 1. The eigenvalue y s = 0.2584 and (S/n) 1 / 2 ~ 1.33 in the core. Note that for a = 1, the 
constraint equation ([31]) for a static core implies that S(0) = 11(0) + ip(0)/Q. So S(0)/II(0) > 1 
for any solution with a static core. The density profile shows a reasonable correspondence with 
the asymptotic behaviour expected taking a = 1; t^(y) oc y . The velocity dispersions and the 
rotation velocity also decrease with decreasing y, but do so a little less rapidly compared to the 
predicted S oc II oc U 2 oc y form. 

In general, we find that, for smaller values of e < 1/6 (or a < 1), while it possible to find 
static core solutions for a sufficiently large S/n ratio, it becomes increasingly difficult to do so 
(obtaining a small enough V m /y m ratio) as one lowers the ratio of S/n, to even slightly smaller 
values. We considered for example a case with e = 1/6, S2 = 2. This turns out to have S/II ~ 1, 
but we found that we could only decrease V m /y m to a value of order unity and get a solution. We 
get V m ~ 1.5 x 10 -2 and dV/dy ~ 0.06 at y m ; so even though the core is not strictly static, the 
LHS of the scaled Euler equation (p0|) is much smaller than each of the individual force terms. 
These nearly cancel each other making the the core quasi-static. We give a plot of all the variables 
for this solution in Figure 5, and a summary of some parameters in Table 1. We see that the shape 
of the density profile for this solution is mid- way between the r~ a oc r~ l form and r~ 2 form. The 
velocity dispersions are decreasing with radius but not as rapidly as predicted for a truly static 
core. 

At this stage it is worthwhile to note the following: Recall that the static core condition 
used to derive the asymptotic scaling properties of the density and velocity dispersions, involves 
assuming not only v r (0) = (the boundary condition adopted above), but also that the radial 
velocity vanishes for a range of radii near the origin. This situation can strictly obtain only if 
particles with a given turn around radius have a minimum radius of approach to the centre; so that 
the core at any radius r is evacuated of particles having turn around radii larger than say, i?t(r). 
Such an "evacuated" core will inturn obtain only if the distribution of angular momentum has a 
"hole" near the origin of {vq^v^) plane. Such an angular momentum distribution is indeed assumed 
(and relevant) in the work of White and Zaritsky (1992). However, in the present work we are 
making the statistical assumption that the distribution of tangential velocities is well described by 
its second moment (viz. the tangential dispersion), thereby excluding distribution functions, which 
have a hole. This assumption is quite reasonable for halo cores with are forming by a general 
inhomogeneous collapse. However, in this case, for any shell of particles which pass the caustic at 
some epoch, there are always some particles with sufficiently small angular momentum, that can 
approach close to the halo core. So the halo core will not be strictly static, a feature which will be 
more and more noticeable, as one decreases tangential velocity dispersions relative to the radial 
dispersions. This may account for our result that (for e < 2/3), as one decreases S/II, the density 
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profile is steeper than the ip oc y a form expected for a strictly static core. 

Finally, for the sake of comparison, we have also looked at numerical examples where 
the tangential velocity dispersions are introduced at the 'turn around' radius (taken to be 
approximately at y = 1) rather than at the shock. In this case one has to solve the moment 
equations numerically, both outside and inside the shock radius, match the solutions across the 
shock, using the shock jump conditions (cf. Eq. ( p5| ) when 111 =0), and find the shock location as 
an eigenvalue to satisfy the boundary conditions in Eq. (^). Figures 6 and 7 give two examples 
with e = 0.4, adopting 11(1) = and E(l) = 0.25 and E(l) = 0.30, respectively. The parameters 
of these solutions are given in Table 1. When E(l) = 0.25, the force due to the tangential velocity 
dispersion at turn around is ~ 13.5% of the radial gravitational force. These examples show very 
similar behaviour to the e = 0.4 solutions discussed above (Figures 2 and 3), where the tangential 
velocities are introduced at the shock. For example, the solution shown in Figure 6, has E/II ~ 1 
in the core. The corresponding density profile is shallower than y~ 2 but steeper than the y~ a 
form, reflecting the fact that only a partial memory of the initial profile is retained by self-similar 
evolution in this case. 

The numerical results of this section shows the importance of tangential velocity dispersions, 
in deciding whether the self similar solution, with an initial density profile shallower than 1/r 2 
(e < 2/3) retains a memory of this initial profile or whether the density profile tends to a universal 
1 /r 2 form. The set of solutions we have given show that for a large enough E /II > 1 , the the core 
density profile is indeed close to the form p oc r~ a , with a = 9e/(l + 3e). For E/II ~ 1, some 
memory of the initial density profile is always retained; the density profile has an asymptotic form 
p oc r~ a , with a < a < 2. When E/II << 1, the density profile goes over to the 1/r 2 form derived 
by FG. Also for shallow initial density profiles with a < 1, one must necessarily have a tangential 
dispersion larger than radial dispersion to get a static core region, retaining the memory of the 
initial density profile. 

5. Discussion and Conclusions 

We have explored here the dynamical restrictions on the structure of dark matter halos 
through a study of cosmological self-similar gravitational collapse solutions, adopting a fluid 
approach to the collisionless dynamics of dark matter. In a companion paper (SC099) we consider 
the possibility that a nested sequence of undigested cores in the center of a halo, which have 
survived the inhomogeneous collapse to form larger and larger objects, determine halo structure 
in the inner regions. For a flat universe with P(k) oc k n , scaling arguments then suggest that the 
core density profile scales as, p oc r~ a with a = a n = (9 + 3n)/(5 + n). However, such arguments 
do not tell us how and in fact whether this form will be realized dynamically. The similarity 
solutions worked out in some detail here, allows us to examine this dynamical issue, in a simple 
tractable manner. 
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The problem of spherical self similar collapse, has often been solved by following particle 
trajectories. We adopted here and in SC099 another approach, examining directly the evolution 
of the moments of the phase space density. For a purely radial collapse, with the initial density 
profile oc r~ 3e , and steeper than r~ 2 , we recover, by demanding that the core be static, the 
asymptotic form of the non-linear density profile: p oc r~ a oc r - 9e /( 1 + 3e ) ( see a l so Padmanabhan 
1996b). For initial density profiles shallower than 1/r 2 , with e < 2/3, we showed that, a static core 
with a non-linear density profile, with a = 9e/(l + 3e), is possible, only if the core has sufficiently 
large tangential velocity dispersions. Infact, one needs Vg > GM/2r. Also if a static core has to 
have a cuspy density profile shallower than 1/r, (with a < 1), one requires Vn > v%- Note that 
when 3e = (3 + n)/2 (as would be relevant for collapse around a typical point in the universe), 
a = a n = (9 + 3n)/(5 + n). 

The consequences of introducing non radial velocity dispersions, in this approach, can only be 
examined in detail, by adopting a closure approximation. In spherical collapse, the skewness of the 
tangential velocities can be assumed to be zero, in the core regions. In fact, in regions where large 
amounts of shell crossing has occurred, one can assume that a quasi "equilibrium" state obtains, 
whereby all odd moments of the distribution function, over (v — v), may be neglected. The radial 
peculiar velocity is then also expected to have negligible skewness, in the core regions. However, 
the radial peculiar velocity will necessarily have a non-zero skewness (non zero third moment) near 
a caustic radius, where collapsing dark matter particles meet the outermost shell of re-expanding 
matter. To take this into account we introduce a fluid approach. In this approach, the effect of 
peculiar velocity skewness is neglected in all regions except at location of the caustic, which we call 
the shock. In the particle picture the shock is where a single stream flow becomes a muti stream 
flow. In the fluid picture it is a where some of the average infall velocity, is converted to velocity 
dispersion. The location of the caustic, y s , in scaled co ordinates, is found as an eigenvalue, to 
the boundary value problem of matching the single stream collapse solution with a core solution, 
adopting V = M = as the boundary condition at y = 0. 

In spherical collapse tangential velocities are only non zero if they are present in the initial 
condition. The shock or the turn around radius, provide a natural location for introducing 
tangential dispersions, into the initial conditions. Our treatment here assumes that the distribution 
of tangential velocities is well described by just its second moment, consistent with the statistical 
assumptions of a quasi-relaxed core. The results of the numerical integration of the moment 
equations, are summarized in Table 1 and are graphically displayed in Figures 1-7. The details of 
one particular solution is also given in Table 2. 

These examples largely bear out the expectations of section 2. First we recover quite well, 
using the fluid approach, the the asymptotic form of the non-linear density profile, for the e = 1 
case, which B85/FG got by solving for the self-similar particle trajectory. Second our solutions 
show the importance of tangential velocity dispersions, in deciding the nature of the core density 
profile, when e < 2/3. In the spherical self similar collapse solutions with e < 2/3, for a large 
enough S/II > 1, one gets p oc r~ a , with a = 9e/(l + 3e). For E/n ~ 1, some memory of the 
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initial density profile is always retained; one gets p oc r~ a , with a < a < 2. When £/n << 1, the 
density profile goes over to the 1/r 2 form derived by FG for radial collapse. Also a < 1, requires 
S/II >> 1, to get a static core region. So if in halo cores tangential velocities are constrained to 
be smaller than radial velocity dispersions, then a cuspy core density profile shallower than 1/r 
cannot obtain, purely by self-similar evolution. 

The results of this work and SC099, illustrate the importance of dynamical considerations 
and hint at features which are likely to obtain in more realistic collapse. If newly collapsing 
material is constrained to mostly contribute to the density at larger and larger radii, then memory 
of initial conditions can be retained. The solutions, with a > 2 (Figure 1), or with a < 2 but a 
large enough tangential dispersion (Figures 2 and 4), illustrate this possibility. However when 
newly collapsing material is able to occupy similar regions as the matter which collapsed earlier, 
the core density profile will only partially reflect a memory of the initial conditions. The solutions 
in Section 4 with a < 2 and S/II ~ 1 (Figures 3, 5 and 6) illustrates this feature. 

In SC099 we have also adopted a complimentary approach, of looking at halo properties 
in numerical simulations of structure formation models having n = —2, —1 and 0. We find that 
the core density profiles of dark matter halos show a large scatter in their properties, but do 
nevertheless appear to reflect a memory of the initial power spectrum (please see SC099 for 
details). The fluid approach adopted here and in SC099 suggests new ways of exploring non linear 
dynamics. Perhaps one can extend analytic approximations like the Zeldovich approximation, 
valid in a single stream flow, to the multi streaming regime, by replacing multistreaming regions 
by regions with velocity dispersions, generated by the Zeldovich type caustics. The fluid approach 
could also be useful to study possible closures of the BBJKY hierarchy. Further one needs to 
extend the self-similar solutions to incorporate a baryonic component; the gas necessarily has an 
isotropic velocity dispersion, and so will have a different dynamical evolution compared to the 
dark matter. We hope to study some of these issues in the future. 

KS thanks Jerry Ostriker and Renyue Cen for an enjoyable collaboration which led to this 
paper. This work was begun when KS visited the Princeton University Observatory, during 
Sept-Nov 1996. Partial travel support to Princeton came from IAU Commission 38. Some of 
the work was done at the University of Sussex where KS was supported by a PPARC Visiting 
Fellowship. He thanks John Barrow, Jerry Ostriker, Ed Turner, the other Princeton and Sussex 
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Fig. 1. — Self similar collapse solution for e = 1 case of B85/FG. The velocity V, scaled density 
tp (solid line in the upper right plot), radial velocity dispersion squared IT, and circular velocity 
squared U 2 are plotted against the scaled radius y. The pre shock spherical infall solution is also 
shown. In the ip-y plot we also show for comparison the density laws ip oc r~ a (dashed line) with 
a = 9e/(l + 3e) and ip oc r~ 2 (dashed -dotted line). 
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Fig. 2. — Self similar collapse solution for e = 0.4, S2 = 0.94. The velocity V, scaled density tp 
(solid line in the upper right plot), radial velocity dispersion squared II (solid line in the lower 
left plot), tangential velocity dispersion squared S (dashed line in the lower left plot) and circular 
velocity squared U 2 are plotted against the scaled radius y. The pre shock spherical infall solution 
is also shown. In the ip-y plot we also show for comparison the density laws ip oc r~ a (dashed line) 
with a = 9e/(l + 3e) and tp oc r~ 2 (dashed -dotted line). 
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Fig. 3. — Self similar collapse solution for e = 0.4, £2 = 0.65. The various quantities shown arc 
same as in Fig. 2. 
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Fig. 4. — Self similar collapse solution for e = 1/6, £2 = 2.2. The various quantities shown are 
same as in Fig. 2. 
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Fig. 5. — Self similar collapse solution for e = 1/6, £2 = 2.0. The various quantities shown are 
same as in Fig. 2. 



-23- 





-2 -1.5 -1 -0.5 ' -2 -1.5 -1 -0.5 

iog(y) log(y) 

Fig. 6. — Self similar collapse solution for e = 0.4, £2(1) = 0.25, 11(1) = 0. Here the tangential 
velocity dispersions have been introduced at the turn around radius corresponding to y = 1. The 
force due to the tangential velocity dispersion at turn around is 13.5% of the radial gravitational 
force. The various quantities shown are same as in Fig. 2. 



-24- 





-2 -1.5 -1 -0.5 -2 -1.5 -1 -0.5 

iog(y) log(y) 

Fig. 7. — Self similar collapse solution for e = 0.4, £2(1) = 0.30, 11(1) = 0. The tangential velocity 
dispersions have been introduced at the turn around radius corresponding to y = 1. The various 
quantities shown are same as in Fig. 2. 
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Table 1. The parameters of the self similar solutions 
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V 2 




t 2 


n 2 


y m 


V m 


i>(y™) 


S(Vm) 


n(wm) 


1 





0.4628 


-0.313 


1.009 





0.433 


1.7E-4 


-4.0E-5 


3.139 





5.464 


0.4 


0.94 


0.4955 


-0.0273 


3.141 


0.94 


0.517 


2.86E-3 


-1.0E-6 


3.641 


1.126 


0.676 


0.4 


0.65 


0.3797 


-0.223 


2.750 


0.65 


0.672 


4.62E-3 


-1.75E-3 


9.691 


4.661 


6.162 


1/6 


2.2 


0.2584 


-0.0473 


6.785 


2.2 


1.232 


1.12E-2 


-2.0E-3 


10.16 


4.021 


2.261 


1/6 


2.0 


0.1999 


-0.163 


7.140 


2.0 


1.584 


1.56E-2 


-1.5E-2 


24.17 


13.45 


9.615 
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0.4 


0.25 


0.2684 


-0.1995 


3.024 


0.78 


0.449 


5.2E-3 


-1.27E-3 


10.67 


4.180 


4.318 


0.4 


0.30 


0.2622 


-0.1521 


3.210 


0.89 


0.363 


4.5E-3 


-4.5E-4 


8.657 


2.909 


2.203 



Table 2. The self similar solutions with e = 0.4, S 2 = 0.94 
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-02 
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.141 





.940 


0.517 


30 


.24 


0.4562 


-2 


.353E- 


•02 


3 


.164 





.948 


0.524 


30 


.39 


0.4226 


-2 


.064E- 


-02 


3 


.184 





.955 


0.530 


30 


,51 


0.3937 


-1 


.832E- 


-02 


3 


.201 





.961 


0.535 


30 


.62 


0.3684 


-1 


.643E- 


-02 


3 


.217 





.966 


0.540 


30 


.72 


0.3265 


-1 


.352E- 


-02 


3 


.243 





.975 


0.548 


30 


.90 


0.2932 


-1 


141E- 


-02 


3 


.264 





.983 


0.555 


31 


.04 


0.2660 


-9 


811E- 


-03 


3 


.283 





.990 


0.560 


31 


.16 


0.2336 


-8 


.036E- 


-03 


3 


.306 





.999 


0.567 


31 


.32 


0.2009 


-6 


.398E- 


-03 


3 


.331 


1 


.008 


0.575 


31 


.49 


0.1710 


-5 


.027E- 


-03 


3 


.355 


1 


.018 


0.583 


31 


.66 


0.1415 


-3 


.798E- 


-03 


3 


.381 


1 


.028 


0.591 


31 


.85 


0.1113 


-2 


671E- 


-03 


3 


.412 


1 


.040 


0.601 


32 


.06 


0.0802 


-1 


.658E- 


-03 


3 


.448 


1 


.054 


0.613 


32 


.32 


0.0506 


-8 


.514E- 


■04 


3 


.492 


1 


.072 


0.626 


32 


.62 


0.0201 


-2 


.262E- 


-04 


3 


.558 


1 


.098 


0.648 


33 


.09 


0.0100 


-8 


.091E- 


-05 


3 


.593 


1 


.113 


0.660 


33 


.33 


0.0080 


-5 


.759E- 


-05 


3 


.603 


1 


.116 


0.663 


33 


.39 


0.0050 


-2 


.493E- 


-05 


3 


.621 


1 


.122 


0.669 


33 


.50 


0.0040 


-1 


.384E- 


-05 


3 


.629 


1 


.124 


0.672 


33 


.54 


0.0032 


-4 


.482E- 


-06 


3 


.637 


1 


.126 


0.675 


33 


.56 


0.0029 


-1 


.000E- 


■06 


3 


.641 


1 


.126 


0.676 


33 
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